
rm(list=ls())

## requires JAGS version 2.1.0, rjags version 2.1.0-2, and R2jags version 0.02-07
## it is incompatible with newer versions of JAGS e.g. 2.2 and 3.x

library(R2jags)
library(foreign)

## dummy variables for missing data
dummies <- TRUE

## choice of violence experience variable
violence.exp <- "individual"

remove.tajiks <- TRUE

load("endorseSetupProTaliban.RData")

indivcov <- cbind(indivcov[,-16], indivcov[,1] * indivcov[,15], indivcov[,3] * indivcov[,15])
colnames(indivcov)[16:17] <- c("pro.taliban.violent.taliban.int", "pro.taliban.violent.ISAF.int")

i <- 1

##########################
##  The Village Model  ##
##########################
data <- list ("N", "J", "K", "Y", "endorser", "village", "district",
              "province", "n.village", "n.district", "n.province", "indivcov",
              "villcov", "distcov")
inits <- function() {list(
                          alpha0 = matrix(seq(-2,2,length.out=(4*J)),J,4), 
                          beta0 = runif(J),
                          x = rnorm(N),
                          s = endorser - 1, 
                          delta.village0 = rnorm(n.village, sd = 5),
                          delta.district0 = rnorm(n.district, sd = 5),
                          delta.cov10 = rnorm(ncol(indivcov), sd = 5),
                          lambda.cov10 = rbind(rep(0, ncol(indivcov)),
                            matrix( rnorm((K-1)*ncol(indivcov), sd = 5),
                                   nrow = K - 1, ncol = ncol(indivcov))),
                          delta.village.cov10 = rnorm(ncol(villcov), sd = 5),
                          delta.district.cov10 = rnorm(ncol(distcov), sd = 5),
                          delta.province = rnorm(n.province, sd = 5),
                          lambda.village0 = matrix(c(rep(0,n.village),
                            rep(0,n.village*(K-1))), n.village, K),
                          theta.village.cov10 = rbind(rep(0, ncol(villcov)),
                            matrix( rnorm((K-1)*ncol(villcov), sd = 5),
                                   nrow = K - 1, ncol = ncol(villcov))),
                          theta.district.cov10 = rbind(rep(0, ncol(distcov)),
                            matrix( rnorm((K-1)*ncol(distcov), sd = 5),
                                   nrow = K - 1, ncol = ncol(distcov))),
                          theta.district0 = matrix(c(rep(0,n.district),
                            rnorm(n.district*(K-1), sd = 5)), n.district, K),
                          theta.province0 = matrix(c(rep(0,n.province),
                            rnorm(n.province*(K-1), sd = 5)), n.province, K),
                          omega = c(0,runif(K-1)),
                          sigma = runif(n.district),
                          psi.district = matrix(c(rep(0,n.district),
                            runif(n.district*(K-1))), n.district, K),
                          psi.province = matrix(c(rep(0,n.province),
                            runif(n.province*(K-1))), n.province, K)
  )}
params <- c("sd.x", "beta", "alpha", "psi.district", "psi.province", "omega", "sigma.district", "sigma.province", "cor.sx", "mean.correct", "delta.village", "delta.district", "delta.province", "lambda.village", "theta.district", "theta.province", "delta.cov1", "lambda.cov1", "delta.village.cov1", "delta.district.cov1", "theta.village.cov1", "theta.district.cov1", "s") 
model <- "ModelVillageDistrictProvinceCovariates.R" 
fit <- jags(data, inits, params, model.file=model, n.iter=5000, n.chains = 1) #takes quite a while -- several days for convergence

sink(paste("results/endorseResultsCovariatesProTal.",violence.exp, "." , i, ".Rout", sep =""))
print(fit)
sink()

save.image(paste("results/endorseResultsCovariatesProTal.",violence.exp, "." , i, ".RData", sep =""))

check <- any(fit$BUGSoutput$summary[, "Rhat"] > 1.2)
counter <- 0
while(check & counter<2000){ 
  fit <- update(fit, n.iter=5000, n.thin=20)
  check <- any(fit$BUGSoutput$summary[, "Rhat"] > 1.2) ## True if any parameters have not reached convergence
  counter <- counter + 1
  
  sink(paste("results/endorseResultsCovariatesProTal.",violence.exp, "." , i, ".Rout", sep =""))
  print(fit)
  sink()
  
  save.image(paste("results/endorseResultsCovariatesProTal.",violence.exp, "." , i, ".RData", sep =""))
  
}  



